人类遗传变异决定24小时节律性基因表达与疾病风险
接收日期：2024年7月23日 作者：Ying Chen(^1), Panpan Liu(^1), Aniko Sabo(^2) & Dongyin Guan (\bigcirc)(^1) 接受日期：2025年4月24日 在线发布日期：2025年5月8日
24小时的生物节律对维持生理稳态至关重要。这些节律的破坏会增加多种疾病的风险。生物节律的形成已知与核心时钟基因的遗传基础有关，但个体遗传变异如何塑造振荡转录组并影响人类时间生理学和疾病风险，目前尚不清楚。本研究通过绘制时间基因表达与基因型的相互作用图谱，确定了影响节律性基因表达的数量性状位点（QTLs）。这些新发现的QTLs被称为节律性QTLs（rhyQTLs），它们揭示了特定基因型人类亚群中未被充分认识的节律性基因。在功能上，rhyQTLs及其相关节律性基因广泛参与重要的时间生理过程，如胆汁酸和脂质代谢。rhyQTLs的发现为基因节律性的遗传机制提供了新见解，揭示了人类疾病风险的变异机制，并为患者的精准时间治疗提供了依据。
24小时的生物节律对维持生理稳态至关重要。这些节律的紊乱会增加患多种疾病的风险。众所周知，生物节律具有由核心时钟基因形成的遗传基础，但个体遗传变异如何塑造振荡的转录转录组，并对人类时间生理和疾病风险做出贡献，在很大程度上是未知的。在这里，我们绘制了时间基因表达和基因型之间的相互作用图谱，以确定与节律性基因表达有关的数量性状位点（qtl）。这些新发现的QTLs被称为节律性QTLs(rhyQTLs)，它们决定了具有特定基因型的人类亚群中以前未被发现的节律性基因。在功能上，rhyQTLs及其相关节律基因广泛参与必要的时间生理过程，包括胆汁酸和脂质代谢。rhyqtl的鉴定揭示了基因节律性的遗传机制，为人类疾病风险的变化提供了机制见解，并为患者提供了精确的时间治疗方法。
生物节律是指周期约为24小时的重复性生理过程。这些节律使哺乳动物能够适应昼夜环境变化（如光暗周期），对维持生理稳态至关重要(^{1,2})。生物节律的破坏已被认为是多种疾病的风险因素，包括2型糖尿病、心血管疾病、消化系统疾病和癌症(^{4-6})。尽管已知不同个体的24小时基因表达和生理过程存在变异(^{7,8})，但这些变异如何影响疾病风险及其机制仍是未解之谜。
此外，目前对节律性基因表达调控机制的理解主要集中在转录因子（TFs）对顺式调控元件（CREs）的作用上，包括核心时钟组件（如BMAL1和REV-ERBs）以及非经典时钟TFs(^{9-11})。动物模型中的遗传或环境干扰可导致或加剧多种疾病(^{12,13})。CREs中的变异与人类基因调控、睡眠障碍、个体时间类型以及其他复杂性状和疾病相关(^{14-16})。尽管这些变异可能影响核心时钟基因的表达，但其具体作用机制尚不明确。此外，遗传变异与特定组织中人类生物节律的关系及其与复杂性状和疾病的联系机制仍未充分探索。由于无法在人类中直接进行遗传操作，本研究利用基因型组织表达（GTEx）项目中的自然遗传变异，分析了其与基因节律性（即基因表达以24小时为周期的周期性波动）的关联，并进一步确定了不同组织中基因节律性与人类表型的关系。
此外，我们目前对节律性基因表达调控机制的理解强调了转录因子（TFs）对顺式调控元件（cre）的调控作用，包括核心时钟元件如BMALl和REV-ERBs，以及非规范时钟TFs-1。在动物模型中，这些调节性tf的遗传或环境扰动可导致或加重多种疾病12,13。假定的cre的变异与人类基因调节、睡眠障碍、个体时间型和其他复杂特征和疾病的变异有关。尽管cre中的这些变异可能影响核心时钟基因的表达，但靶向这些cre的特异性tf及其作用机制仍不清楚。此外，遗传变异与特定组织中人类生物节律之间的关系，以及这些变异与人类复杂性状和疾病之间的联系机制，在很大程度上尚未得到探索。由于直接遗传操作在人类中是不切实际的，我们在基因型组织表达（GTEx）项目中研究自然遗传变异，以确定基因节律性扰动的关联，基因节律性是指基因表达以24小时为周期的周期性波动，并进一步确定基因节律性与各种组织中人类表型之间的关系。
Genotype Tissue Expression (GTEx) Project
结果
我们首先利用GTEx数据在838名个体的45种组织中建立了约3900万对遗传变异-基因关联（补充数据1 ）(^{17})，然后评估了遗传变异与其附近基因节律性表达的关联（图1a）。通过谐波回归分析每种遗传变异-基因对在特定基因型亚群中的节律性，我们发现，45种组织中共有中位数为3200个基因在至少一个基因型亚群中呈现节律性表达。
补充图1：多个测试更正和跨组织的rhyQTL发现。a每个基因在45个组织中独立测试的估计有效数量（Mef）的分布。计算Mer以确定有效独立测试的数量，为估计I型错误（假阳性）和应用多个测试更正提供基础。b从表型洗牌数据中发现的rhy变异和rhyqtl的数量。蓝色虚线表示在真实数据集中识别的rhy变体的5%。红色虚线表示在真实数据集中发现的rhyQTLs的10%。c在45个组织中发现的rhyQTLs的数量。rhyQTL编号是根据独立LD块的数量来定义的，其中遗传变异被分组到不同的LD块中。
例如，心脏组织中与心脏移植发展相关的同种异体炎症因子1（AIF1）的表达仅在SNP rs7740525的GG基因型亚群中呈现节律性（图1b）。值得注意的是，AIF1的整体表达水平在这三种基因型中无显著差异。这表明SNP rs7740525与AIF1节律性的关系无法通过表达数量性状位点（eQTLs）解释 ，后者与不同基因型间基因表达量的变异相关。我们将与24小时节律性基因表达变异相关的遗传位点称为节律性数量性状位点（rhyQTLs），将至少与一个rhyQTL相关的基因称为节律性基因（rhyGenes） 。
We first established a median of 39 million genetic variant-gene pairs across 45 tissues from 838 individuals using GTEx data (Supplementary Data 1)17 and then assessed the association between genetic variants and rhythmic expression of their nearby genes (Fig. 1a). Using harmonic regression to evaluate the gene rhythmicity in a subpopulation with specific genotype for each genetic variant-gene pair, we identified a median of 3200 genes across 45 tissues that are rhythmically expressed in at least one genotype subpopulation. A
图1 | 人类组织中节律数量性状位点（rhyQTLs）的全基因组定位 a. rhyQTL定位和rhyGene鉴定的研究设计 该研究分析了45种人体实体组织 ，通过多步骤流程鉴定影响基因表达节律性的遗传变异（rhyQTLs）和节律性表达基因（rhyGenes）。图中显示的数字（基因数量、基因-变异位点对）是这45个组织的中位值。
关键步骤解析
步骤3 ：使用**谐波回归（harmonic regression）**评估两个关键参数：
p值 （统计显著性）振幅（amplitude） （峰值与谷值的log2倍数变化）
步骤4 ：通过**谐波方差分析（harmonic ANOVA, HANOVA）**评估基因节律性的变异。
进一步使用Benjamini-Hochberg (BH)校正 计算调整后的p值，控制假阳性率。
示例说明
步骤三 ：发现基因 AIF1 在GG基因型群体中表达呈现昼夜节律（p=4.8×10⁻⁴，振幅=0.92），但在GT/TT群体中无节律。步骤四 ：通过HANOVA证明 AIF1 的节律性差异是由SNP rs7740525的基因型（GG vs. GT/TT）导致的，因此rs7740525被定义为 AIF1 的 rhyQTL 。
简言之，步骤三筛选节律基因，步骤四鉴定调控这些节律的遗传变异 。
可视化工具 该图使用 BioRender 绘制（作者：Guan, D. (2025)，链接：BioRender.com/nlykndj ）。
核心概念
rhyQTL（节律数量性状位点） ：影响基因表达昼夜节律的遗传变异（如SNPs）。rhyGene（节律性基因） ：呈现周期性（如昼夜节律）表达模式的基因。谐波回归 & HANOVA ：用于分析周期性数据的统计方法，比传统ANOVA更适合节律研究。
该研究为理解基因表达节律的遗传调控 提供了全基因组层面的视角。
补充图2a
补充图2b-c 补充图3
由于特定基因型亚群中的节律性表达模式可能被其他基因型亚群中的非节律性表达掩盖，某些基因的节律性可能仅在特定基因型亚群中显现，而在整体人群中无法观察到。事实上，我们的rhyQTL分析揭示的节律性基因数量是先前未考虑遗传变异的全人群研究中发现的4倍（图1c ）。例如，在心脏组织中，新发现了2748个节律性基因，其节律性表达仅存在于特定基因型中 。即使在整体人群中检测到的681个节律性基因中，98%受rhyQTLs影响，仅2%的节律性基因独立于附近遗传变异的顺式调控 （图1d）。图1e ）。13个不受遗传变异影响的节律性基因包括对心脏功能至关重要的基因，如心肌细胞特异性表达的肌球蛋白轻链激酶3（MYLK3），其节律性可能是进化保守的。
在不同组织中，rhyGenes及其相关变异（rhyVariants）的数量存在差异，例如内脏脂肪组织中有3965个rhyGenes和183,758个rhyVariants，而黑质中仅有202个rhyGenes和653个rhyVariants（图2a及补充数据2 ）。补充图1c ）。通过随机抽样和昼夜RNA-seq数据验证，我们确认了GTEx数据在检测人类节律性表达变异相关遗传变异方面的有效性（补充图4-5 ）。
除顺式调控外，我们还探索了rhyQTLs对节律性基因表达的远程（trans）效应。以核心时钟基因为例，我们发现其节律性表达与多种组织中的遗传变异相关（补充图6a ）。
通过量化rhyGenes在组织中的分布，我们发现了793个在20种以上组织中常见的rhyGenes，这些基因与褪黑素、失眠和睡眠调节等普遍昼夜功能相关（图2b ）。
进一步分析发现，仅3-37%的rhyQTLs同时也是eQTLs（图3a），表明rhyQTLs主要调控基因节律性而非整体表达水平。基因组功能注释显示，rhyQTLs在增强子区域的富集程度高于eQTLs（图3b及补充数据5），这与时钟调节因子偏好结合增强子的发现一致(^{26,29})。通过扫描rhyQTLs中的人类TF结合基序，我们发现脑组织中rhyQTLs高度富集核心时钟基因（如BMAL1、CLOCK和NR1D1）的结合基序 （图3c-d及补充数据6），表明遗传变异与时钟调节因子在调控生物节律变异中的相互作用。
为解析rhyQTLs对人类复杂性状和疾病变异的贡献，我们从GWAS目录中检索了约50万SNP-性状/疾病关联，评估了不同组织中rhyQTLs在GWAS标签SNPs中的富集情况（补充图10a-c及补充数据7 ）。结果显示，肝脏rhyQTLs在GWAS标签SNPs中富集程度最高（图4a），尤其与脂质/脂蛋白标志物和肝酶水平等肝脏相关功能直接相关的性状（图4b）。通过分层连锁不平衡评分回归 （LDSC）分析15种脂质或脂蛋白相关性状，发现rhyQTLs和eQTLs分别解释了中位数15.8%和22.9%的SNP遗传力（图4c）。例如，LGR4的rhyQTLs与亚油酸与总脂肪酸比例的GWAS信号共定位，而其eQTLs无此关联（图4d）。类似地，EHMT2和SELENON的rhyQTLs分别与低密度脂蛋白胆固醇和总胆固醇水平的GWAS信号共定位（图4e-f）。
在肥胖易感人群的GWAS分析中，我们发现了186个SNP-性状关联，包括染色体11上TNFα水平的GWAS信号（图4g及补充数据10-11）。这些SNPs在多种组织（尤其是肝脏）中作为rhyQTLs调控IGF2的节律性表达，而IGF2通过调节TNFα信号介导炎症反应(^{48,49})。富集分析表明，肝脏、结肠和脂肪组织中的rhyQTLs对肥胖相关性状有显著贡献（图4h），提示节律性遗传调控在肥胖相关代谢功能障碍中的潜在作用。
讨论
通过分析全基因组基因型和附近基因的节律性，我们绘制了一种新型QTL——rhyQTL，它决定了多种人类组织中节律性基因表达的变异。rhyQTLs为个体时间生理学变异（如胆汁酸和胆固醇代谢）提供了新见解。在基因组坐标和调控机制上，rhyQTLs与eQTLs存在显著差异。rhyQTLs更富集于增强子区域，这与基因节律性表达依赖于环境线索和内在时钟机器的概念一致(^{9,20})。功能上，rhyQTLs（独立于eQTLs）广泛贡献于人类疾病风险变异，如代谢紊乱和心血管疾病。总之，rhyQTLs介导的节律性基因表达变异为GWAS信号提供了新解释，为理解人类时间生理学的遗传变异提供了新视角，并为基于患者遗传背景优化时间治疗提供了依据。
尽管本研究具有重要价值，但仍存在一些局限性。与传统QTL分析不同，rhyQTL分析缺乏效应大小估计（如beta值或Z分数），限制了精细定位工具的应用。未来需开发新统计方法以改进因果变异的识别。此外，本研究 仅在内脏脂肪组织中探索了rhyQTLs的性别特异性，未来需在更多组织中扩展分析 。
方法
步骤一
数据获取与预处理（用于rhyQTL定位，图1a步骤1）
基因型数据获取与质量控制
数据来源：从dbGaP数据库（编号phs000424.v8.p2）下载GTEx项目中 838名个体 的基因型数据。（SNP分型、各组织表达量、不同时间的表达量）
tabix 支持多种常见的基因组学文件格式，包括：
VCF （Variant Call Format）：用于存储基因组变异信息。BED （Browser Extensible Data）：用于描述基因组特征的区间。GFF/GTF （General Feature Format / Gene Transfer Format）：用于描述基因组注释信息。
质控步骤：
使用全基因组关联分析工具 PLINK 过滤遗传变异，保留以下变异：（standard quality control criteria53.）
最小等位基因频率（MAF）≥ 0.01（排除罕见变异）。
哈迪-温伯格平衡（HWE）p值 ≥ 10⁻⁶（排除偏离群体遗传平衡的变异）。
仅保留常染色体变异（排除性染色体和线粒体变异）。
保留以下组织：tissues with more than 100 samples were retained, as indicated in Supplementary Data 1.仅保留 样本量超过100例 的组织
QTL analysis协变量 。这些协变量通过线性回归模型进行回归分析，以最小化混杂效应并提高分析准确性
考虑到batch effects and population structure in rhyQTL mapping across tissues。5个基因型主成分（PCs）、测序平台、文库构建协议和供体性别（遵循GTEx联盟标准流程）。
five genotype principal components (PCs),
sequencing platform,
library construction protocol, and
donor sex17.
此外，我们还纳入了缺血时间、年龄和死亡类型等参数，这些数据被用于在全人群范围内识别节律基因，而不考虑基因型效应。（ 额外校正 缺血时间 （样本采集后至冷冻保存的时间）、供体年龄 和 死亡原因 ，以控制其对基因表达节律的潜在影响。 ）21.本研究中的时间定义为供体内在昼夜节律相位，而非死后样本采集时间。这一内在相位整合了来自同一捐赠者的组织样本的相对昼夜节律相位，其基础假设是这些相位在不同捐赠者间具有保守性，且不受确切死亡时间的影响17。（重要的是，死亡时间或样本冷冻保存时间本身都不足以界定昼夜节律性，因为捐赠者的时间型、所处时区位置以及死亡原因等因素都可能掩盖真实的昼夜节律相位。） 仅依赖死亡时间或冷冻时间无法准确反映真实的昼夜节律（需整合供体的时间类型等复杂因素）。
关键概念 ：“时间” 定义为 供体内在的昼夜节律相位 （而非死亡时间或样本采集时间）。原理 ：GTEx数据集中收集了54个组织的RNA-seg数据。本研究排除了6个少于100个基因型样本的组织：肾髓质、输卵管、宫颈外、宫颈内、膀胱和肾皮质。此外，还排除了培养成纤维细胞ebv转化淋巴细胞系和全血。结果，共纳入45个实体组织进行rhyQTLs的鉴定。对于每个剩余的组织，只纳入以下标准的样本：RNA完整性数（RIN）大于6，总reads数超过4000万，唯一映射的reads比例大于0.8，阳性缺血时间和自溶评分为1或0。此外，每个组织中样本中平均计数小于10的基因被排除在外
注意事项
数据访问权限 ：时间信息处理 ：内在节律相位的推断需依赖额外数据（如供体的活动记录或生物标志物），文中未公开具体算法。 代码复现 ：rhyQTL_mapping.sh）和预处理脚本。
此部分为rhyQTL分析的基础步骤，后续需进行遗传变异-基因对的节律性建模（谐波回归）和差异检验（HANOVA）。
英文术语
中文翻译
说明
Minor Allele Frequency (MAF) 最小等位基因频率
群体中次常见等位基因的出现频率，MAF < 0.01的变异通常视为罕见变异。
Hardy-Weinberg Equilibrium (HWE) 哈迪-温伯格平衡
群体遗传学中基因型频率的稳定状态，偏离可能提示技术误差或自然选择。
Principal Components (PCs) 主成分
用于校正群体分层（如种族差异）的统计变量。
Ischemic time 缺血时间
样本从死亡到冷冻保存的时间间隔，可能影响RNA降解速率。
Internal circadian phase 内在昼夜节律相位
基于供体生理节律（如睡眠-觉醒周期）推断的时间点，而非实际采样时间。
步骤二
遗传变异-基因配对及节律性评估（图1a中的步骤2和3） 步骤2：建立遗传变异-基因对
顺式调控变异筛选 ：
对于每个蛋白质编码基因，选择其转录起始位点（TSS）上下游 ±1 Mb 范围内的遗传变异 （如SNP）。 SNP是最常见的顺式变异形式 。
该范围是 顺式QTL（cis-QTL） 研究的常用标准。
统计结果 ：
在45个组织中，平均每个基因对应 2.9k个顺式遗传变异。
总共建立了约 3900万（39 million）个变异-基因对 。
数据过滤 ：
保留至少在 2或3种不同基因型亚群 中样本量 > 50可分析的变异-基因对，以确保统计稳健性。
步骤3：评估基因表达的节律性
节律性分析方法 ：
采用 谐波回归（cosinor回归） ，拟合基因表达的昼夜波动模式。
模型公式 ：
E x y = μ x y + a x y cos ( 2 π T t y ) + b x y sin ( 2 π T t y ) + ϵ
(E_{xy})：基因 (x) 在个体 (y) 中的表达水平。
(\mu_{xy})：基因 (x) 在基因型 (y) 中的平均表达量。
(a_{xy}, b_{xy})：余弦和正弦项的系数。
(T = 24) 小时（昼夜节律周期）。
(\epsilon)：随机误差（噪声）。
计算振幅（amp）和相位（phase） ：
amp = 2 a x y 2 + b x y 2
Misplaced & 3. **统计检验**： - 使用 **似然比检验（LRT）** 比较 **谐波模型（有节律）** vs. **平坦模型（无节律，即 \(a_{xy} = b_{xy} = 0\)）**。 - 在文章中提到的似然比检验（Likelihood Ratio Test, LRT）中，对数似然值的计算是通过比较两个模型的对数似然值来完成的：一个模型是包含所有参数的完整模型（Alternative Model），另一个是简化模型（Null Model），其中一些参数被固定为特定值（通常是0）。 对于给定的基因表达数据，我们使用谐波回归（Harmonic Regression）来评估基因表达的节律性。谐波回归模型可以表示为： \[ E_{xy} = \mu_{xy} + a_{xy} \cos\left(\frac{2\pi t_y}{T}\right) + b_{xy} \sin\left(\frac{2\pi t_y}{T}\right) + \epsilon \] 其中，\(E_{xy}\) 是基因 \(x\) 在个体 \(y\) 中的表达水平，\(\mu_{xy}\) 是基因 \(x\) 在个体 \(y\) 中的平均表达水平，\(a_{xy}\) 和 \(b_{xy}\) 是余弦和正弦项的系数，\(t_y\) 是个体 \(y\) 的时间点，\(T\) 是周期（在本研究中为24小时），\(\epsilon\) 是随机误差。 对数似然值（Log-likelihood）是通过最大化似然函数来计算的，似然函数是给定参数下观测数据的概率。在谐波回归中，对数似然值可以通过最小化残差平方和来计算，因为假设误差项 \(\epsilon\) 服从正态分布。 对于完整模型（包含所有参数的模型），对数似然值为： \[ \ell_{\text{full}} = -\frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (E_{xy} - \hat{E}_{xy})^2 \] 其中，\(n\) 是观测数，\(\sigma^2\) 是误差项的方差，\(\hat{E}_{xy}\) 是模型预测的基因表达水平。 对于简化模型（参数被固定的模型），对数似然值为： \[ \ell_{\text{null}} = -\frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (E_{xy} - \bar{E}_{xy})^2 \] 其中，\(\bar{E}_{xy}\) 是简化模型预测的基因表达水平，通常是基因表达的平均值。 似然比检验统计量（Likelihood Ratio Statistic）是两个对数似然值的差： \[ \text{LRT} = -2 (\ell_{\text{null}} - \ell_{\text{full}}) \] 在零假设下，这个统计量服从自由度为参数差的卡方分布。通过比较这个统计量和卡方分布的临界值，我们可以确定是否拒绝零假设，即是否存在基因表达的节律性。 在实际应用中，对数似然值的计算通常由统计软件或编程语言（如R或Python）中的回归分析函数自动完成。在本研究中，使用了R语言中的`lrtest`函数来计算似然比检验统计量。 - 计算 **p值**（R语言 `lrtest` 函数），评估==节律性是否显著==。 4. **多重检验校正**： - 由于每个基因对应多个变异，需校正假阳性。（即虽有节律性，但并非某一变异带来，而是该基因对应的另一个变异） - 考虑到 **连锁不平衡（LD）**，计算每个基因的 **有效独立变异数（Meff）**。 . Adjusting p values based on Meff可以校正假阳性 - 校正阈值： - 中位数 **Meff = 15**（每个基因）。 0.01/15 约为5 × 10-4，所以 - 设定 **显著性阈值 \(p \leq 5 \times 10^{-4}\)**（相当于Bonferroni校正 \(0.01/15\)）。 - **验证该阈值合理性**： - 对基因型数据进行随机置换（shuffle）随机打乱，发现假阳性率 **< 5%**，支持该阈值的可靠性。 5. **最终筛选标准**： - 保留满足以下条件的变异-基因对（**潜在节律变异-rhyVariant & 节律基因-rhyGene**）： - \(p \leq 5 \times 10^{-4}\) - **峰值-谷值倍数变化（fold change）≥ 1.5**（至少一种基因型）。 --- ### **关键点总结** - **顺式变异范围**：TSS ±1 Mb（涵盖增强子、启动子等调控元件）。 - **节律性分析**：基于 **24小时周期** 的谐波回归，计算振幅和相位。 - **统计严格性**：通过 **Meff校正** 和 **置换检验** 控制假阳性。 - **最终输出**：显著关联的 **rhyQTLs**（影响基因昼夜表达的遗传变异）。 ## 步骤四 节律性变异的评估（图1a步骤4）： 使用dryR包中的drylm函数和贝叶斯信息准则（BIC）评估不同基因型间节律性差异。通过20次降采样和G检验确定差异节律性的显著性，最终定义==rhyVariants的标准==包括：（1）至少在一个基因型亚群中p≤5×10\(^{-4}\)且振幅≥1.5；（2）dryR和G检验显示显著差异节律性；（3）校正后q(BH)<0.05。通过LD区块划分确定独立rhyQTLs。 性能评估：通过随机抽样验证rhyQTLs检测方法的稳健性（补充图4），并利用昼夜RNA-seq数据验证GTEx中鉴定的rhyQTLs（补充图5）。此外，通过NR1D1相关rhyQTLs分析了其对远端基因节律性表达的远程效应（补充图6）。 通路富集分析：使用enrichR包和Enrichr数据库进行，背景基因列表为相应组织中表达的基因（平均计数>10）。显著性阈值设为校正后p<0.05。 GWAS富集分析：从GWAS目录获取SNP-性状/疾病关联，过滤后保留242,822个高质量关联。通过随机匹配SNP数量和MAF分布计算rhyQTLs在GWAS标签SNPs中的富集程度。 LDSC回归：分析15种脂质或脂蛋白相关性状中rhyQTLs和eQTLs的遗传力贡献比例。 共定位分析：使用COLOC包评估rhyQTLs与GWAS信号的共定位，后验概率PP4>0.70视为显著。 VIVA研究GWAS：对815名肥胖易感儿童进行75种肥胖相关性状的GWAS分析，鉴定显著关联（p≤7.7×10\(^{-8}\)）并评估rhyQTLs的富集。 #### 数据与代码可用性 GTEx项目V8的原始基因组测序数据可通过dbGaP获取（phs000424.v8.p2）， 基因表达数据从GTEx门户网站下载。 rhyQTLs的汇总统计信息见补充数据2， 自定义源代码发布于GitHub（https://github.com/YingChen10/rhyQTL）\(^{82}\)。 您的理解基本正确，但需要根据研究目标决定是否需要下载两个网站的数据。以下是具体建议： ### **1. 数据需求分析** | 数据类型 | 存放位置 | 是否需要下载？ | |------------------------|------------------------------|------------------------------------------------------------------------------| | **原始基因组测序数据** | dbGaP (phs000424.v8.p2) | 仅当您需要：<br>• 重新分析原始序列（如call变异）<br>• 验证特定样本的测序质量 | | **基因表达矩阵** | GTEx Portal (bulk_tissue_expression) | 如果只需分析基因表达量、eQTL/rhyQTL等，**只需下载此数据**（已处理好的标准化表达量） | ### **2. 注意事项** #### **（1）GTEx Portal表达数据的优势** - **已预处理**：包含TPM/FPKM标准化表达矩阵，可直接用于分析。 - **样本丰富**：涵盖838名捐赠者的45种组织，附带临床元数据（年龄、性别、死亡时间等）。 - **免审批**：表达数据为开放访问（Open Access），无需申请权限。 #### **（2）dbGaP原始数据的适用场景** - **受限访问**：需提交研究提案并通过dbGaP审批（涉及人类数据隐私）。 - **使用场景举例**：<br> • 开发新的RNA-seq分析流程<br> • 探索GTEx未提供的转录本异构体或罕见变异<br> • 整合其他组学数据（如外显子组测序） --- ### **3. 推荐下载策略** #### **如果您的目标是：** - **分析基因节律性（如本文）** → **仅需GTEx表达数据**（含时间信息文件 [Zenodo链接](https://doi.org/10.5281/zenodo.7215362)）。 - **挖掘新的遗传变异** → 需额外下载dbGaP原始数据（需审批）。 #### **GTEx表达数据下载步骤：** 1. 访问 [GTEx下载页面](https://www.gtexportal.org/home/downloads)。 2. 选择 **"Bulk Tissue Expression"** → 下载 **"Gene TPMs"**（推荐）或 **"Gene Counts"**。 3. 同步下载 **样本属性表**（含组织类型、供体信息等）。 --- ### **4. 数据关联示例** 若您研究rs7740525对AIF1基因的调控： - **GTEx数据**：提供AIF1在各组织中的表达量，验证其是否在心脏中呈现节律性。 - **dbGaP数据**：如需分析该SNP附近的测序覆盖度或寻找共定位变异，则需原始BAM文件。 --- ### **总结** - **常规表达分析**：GTEx Portal数据足够（更快捷）。 - **深度遗传分析**：需结合dbGaP原始数据（流程复杂）。 - 建议先尝试用GTEx数据完成初步分析，必要时再申请dbGaP权限。 #### 参考文献（略） #### 致谢 感谢所有提供宝贵意见和技术支持的同事，以及资助本研究的相关机构。 #### 作者贡献 Y.C.和D.G.设计研究、解读数据并撰写论文，所有作者修订并批准最终稿。Y.C.和P.L.负责基于RNA-seq数据的SNP calling，Y.C.和A.S.处理VIVA项目数据。 #### 竞争利益 作者声明无竞争利益。 #### 附加信息 补充信息及通讯材料请求详见原文链接。 文件名：补充数据1-11 补充数据1. 本研究用于rhyQTL定位的组织样本汇总，包括组织名称及样本量。 补充数据2. 45种组织中rhyQTL的统计摘要。该数据已上传至Zenodo平台（DOI: 10.5281/zenodo.15127913），包含多个按"组织名称.txt"命名的文本文件。每份文件包含19列数据，具体如下： 1. ID - rhy变异体标识符 2. Chromosome - rhyQTL所在染色体 3. Position - rhyQTL位点坐标 4. REF - 参考等位基因 5. ALT - 替代等位基因 6. rhyGene.ID - 节律基因ID 7. rhyGene.name - 节律基因名称 8-10. Sample.size.0/1/2 - 按样本量排序的三个基因型亚组样本量 8. pval_0 - 第一基因型亚组的p值（11-16列为谐波回归拟合参数） 9. phase_0 - 第一亚组的表达峰值相位（24小时制） 10. amp_0 - 第一亚组的对数化峰谷振幅值（log2） 11. pval_1 - 第二基因型亚组的p值 12. phase_1 - 第二亚组的表达峰值相位 13. amp_1 - 第二亚组的对数化峰谷振幅值 14. Chosen.model - DryR选择的模型（取值为2/3/5的整数） 15. G.test_pval - G检验p值 16. pval - 基因型组间节律表达差异的统计学显著性p值 补充数据3. 脑组织中NR1D1相关顺式rhyQTL介导的反式节律基因。 补充数据4. 节律基因峰值相位。每个节律基因的峰值相位定义为表达达到峰值的时间点。若基因存在多个rhyQTL，则取所有对应相位的中位数作为该基因相位。 补充数据5. rhyQTL在基因组调控元件中的富集分析。通过计算观察值与预期值的比值比（odds ratio）评估富集程度，基线区域为人类基因组中所有功能及推定功能区域的联合集。 补充数据6. 45种组织中rhyQTL位点的转录因子结合基序富集分析。通过构建2×2列联表分析转录因子结合基序与遗传变异的共现关系，使用Fisher精确检验计算比值比和p值评估富集显著性。文件按"组织名称.txt"格式存储，包含9列数据： 1. Motif - 基序序列 2. Name - 最佳匹配基序及详细信息 3. TF - 转录因子名称 4-7. 分别表示"rhyQTL+/基序+"、"rhyQTL+/基序-"、"rhyQTL-/基序+"、"rhyQTL-/基序-"的SNP数量 4. OddsRatio - Fisher精确检验计算的比值比 5. p.value - 富集分析显著性p值 补充数据7. rhyQTL相关性状/疾病。列出GWAS中鉴定到≥3个与rhyQTL重叠SNP的性状/疾病。 补充数据8. 人类性状/疾病分类摘要。GWAS Catalog中的性状/疾病按其生理功能分为13大类。 补充数据9. 与15种脂质/脂蛋白相关性状共定位的rhyQTL节律基因。使用COLOC PP4统计方法鉴定与GWAS信号共享变异位点的rhyQTL，并列出影响特定人类表型的节律基因。 补充数据10. VIVA队列中测量的肥胖相关性状。 补充数据11. GWAS鉴定的肥胖相关性状SNP。基于75个肥胖相关表型，按照标准GWAS分析流程（https://github.com/MareesAT/GWA_tutorial）进行分析。 （注：专业术语处理说明： 1. rhyQTL/rhyGene等专业缩写保留不译，通过上下文明确含义 2. DryR/COLOC等分析方法名称保留原名 3. 统计学术语如odds ratio统一译为"比值比" 4. 技术术语如"peak phase"译为"峰值相位"以符合生物学表述习惯 5. 长难句采用拆分法处理，如补充数据6的描述拆分为流程说明和文件结构两部分） File Name: Supplementary Data 1-11 Supplementary Data 1. Summary of tissues used for rhyQTL mapping in this study, including their names and sample sizes. Supplementary Data 2. Summary statistics of rhyQTLs in 45 tissues. The Supplementary Data 2 has been uploaded to Zenodo: https://doi.org/10.5281/zenodo.15127913. The folder contains multiple .txt files named following the pattern tissue_name.txt. Each file includes data organized into 19 columns, which are listed as follows: 1. ID - Identifies the rhyVariant. 2. Chromosome - Indicates the chromosome of rhyQTL. 3. Position - Specifies the position of rhyQTL. 4. REF - Reference allele of rhyQTL. 5. ALT - Alternative allele of rhyQTL. 6. rhyGene.ID - Gene id for the rhyGene. 7. rhyGene.name - Name of the rhyGene. 8-10. Sample.size.0/1/2 - Sample size of the three genotype subgroups, ordered by sample size. 11. pval_0 - p value for the first genotype subgroup. Columns 11-16 are parameters from the harmonic regression fit. 12. phase_0 - Peak phase (24h scale) for the first genotype subgroup. 13. amp_0 - log2 amplitude peak-to-trough for the first genotype subgroup. 14. pval_1 - p value for the second genotype subgroup. 15. phase_1 - Peak phase (24h scale) for the second genotype subgroup. 16. amp_1 - log2 amplitude peak-to trough for the second genotype subgroup. 17. Chosen.model - The selected model given by DryR, with an integer value of 2, 3, or 5. 18. G.test_pval - p value from the G-test. 19. pval - p value reflects the statistical significance of differential rhythmic expression between the two genotype groups. Supplementary Data 3. Trans-rhyGenes mediated by NR1D1-associated cis-rhyQTLs in brain tissues. Supplementary Data 4. Peak phases of rhyGenes. The peak phase of each rhyGene is defined as the time point at which it exhibits the expression peak. For rhyGenes with more than one rhyQTLs, the median value of phases corresponding to all rhyQTLs is taken as the phase of the rhyGene. Supplementary Data 5. Enrichment of rhyQTLs in genomic regulatory elements. The enrichment is calculated as an odds ratio score using the numbers of observed QTLs and expected QTLs located in each annotation category compared to those in baseline regions. Baseline regions are unions of all functional and putative functional regions within the human genome. Supplementary Data 6. TF motif enrichment in genomic loci harboring rhyQTLs across 45 tissues. The relationship between the cooccurrence of transcription factor (TF) binding motifs and genetic variants at different genomic coordinates was analyzed. Genetic variants were categorized based on whether they were identified as rhyQTLs and whether they were located within the motif, forming a 2 × 2 contingency table for each motif. The folder contains multiple .txt files named following the pattern tissue_name.txt. Each file includes data organized into 9 columns, which are listed as follows: 1. Motif: The specific motif sequence being analyzed. 2. Name: The best matching motif and additional details about the motif 3. TF : The name of the transcription factor associated with the motif. 4. rhyQTL+Motif+: The number of SNPs that are both rhyQTLs and located within a motif binding site. 5. rhyQTL+Motif-: The number of SNPs that are rhyQTLs but are not located within a motif binding site. 6. rhyQTL-Motif+: The number of SNPs that are not rhyQTLs but are located within a motif binding site. 7. rhyQTL-Motif-: The number of SNPs that are neither rhyQTLs nor located within a motif binding site. 8. OddsRatio: The odds ratio calculated using Fisher's exact test, indicating the strength of the enrichment. 9. p.value: The p value calculated using Fisher's exact test, indicating the statistical significance of the enrichmentodds ratio was calculated to assess rhyQTL enrichment within motifs, and statistical significance was determined using a two-sided Fisher's exact test. Supplementary Data 7. rhyQTL-associated traits/diseases. Traits/diseases with more than three associated SNPs identified in GWAS that also identified as rhyQTLs are listed. Supplementary Data 8. Summary of human trait or disease categories. Traits/diseases from GWAS Catalog are classified into 13 categories based on their specific physiological functions. Supplementary Data 9. rhyGenes with rhyQTLs that colocalize with loci linked to 15 lipid/lipoprotein-related traits. rhyQTLs sharing variants with significant GWAS signals in 15 lipid/lipoprotein-related traits were identified using COLOC PP4 statistics to estimate their effects on specific human phenotypes. rhyGenes whose rhyQTLs colocalized with GWAS signals for these traits are listed. Supplementary Data 10. Obesity-related traits measured in VIVA cohort. Supplementary Data 11. SNPs associated with obesity-related traits identified by GWAS. GWAS analysis on these 75 obesity-related phenotypes was conducted following the GWAS standard pipeline https://github.com/MareesAT/GWA_tutorial. 3. **统计检验**： - 使用 **似然比检验（LRT）** 比较 **谐波模型（有节律）** vs. **平坦模型（无节律，即 \(a_{xy} = b_{xy} = 0\)）**。 - 在文章中提到的似然比检验（Likelihood Ratio Test, LRT）中，对数似然值的计算是通过比较两个模型的对数似然值来完成的：一个模型是包含所有参数的完整模型（Alternative Model），另一个是简化模型（Null Model），其中一些参数被固定为特定值（通常是0）。 对于给定的基因表达数据，我们使用谐波回归（Harmonic Regression）来评估基因表达的节律性。谐波回归模型可以表示为： \[ E_{xy} = \mu_{xy} + a_{xy} \cos\left(\frac{2\pi t_y}{T}\right) + b_{xy} \sin\left(\frac{2\pi t_y}{T}\right) + \epsilon \] 其中，\(E_{xy}\) 是基因 \(x\) 在个体 \(y\) 中的表达水平，\(\mu_{xy}\) 是基因 \(x\) 在个体 \(y\) 中的平均表达水平，\(a_{xy}\) 和 \(b_{xy}\) 是余弦和正弦项的系数，\(t_y\) 是个体 \(y\) 的时间点，\(T\) 是周期（在本研究中为24小时），\(\epsilon\) 是随机误差。 对数似然值（Log-likelihood）是通过最大化似然函数来计算的，似然函数是给定参数下观测数据的概率。在谐波回归中，对数似然值可以通过最小化残差平方和来计算，因为假设误差项 \(\epsilon\) 服从正态分布。 对于完整模型（包含所有参数的模型），对数似然值为： \[ \ell_{\text{full}} = -\frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (E_{xy} - \hat{E}_{xy})^2 \] 其中，\(n\) 是观测数，\(\sigma^2\) 是误差项的方差，\(\hat{E}_{xy}\) 是模型预测的基因表达水平。 对于简化模型（参数被固定的模型），对数似然值为： \[ \ell_{\text{null}} = -\frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (E_{xy} - \bar{E}_{xy})^2 \] 其中，\(\bar{E}_{xy}\) 是简化模型预测的基因表达水平，通常是基因表达的平均值。 似然比检验统计量（Likelihood Ratio Statistic）是两个对数似然值的差： \[ \text{LRT} = -2 (\ell_{\text{null}} - \ell_{\text{full}}) \] 在零假设下，这个统计量服从自由度为参数差的卡方分布。通过比较这个统计量和卡方分布的临界值，我们可以确定是否拒绝零假设，即是否存在基因表达的节律性。 在实际应用中，对数似然值的计算通常由统计软件或编程语言（如R或Python）中的回归分析函数自动完成。在本研究中，使用了R语言中的`lrtest`函数来计算似然比检验统计量。 - 计算 **p值**（R语言 `lrtest` 函数），评估==节律性是否显著==。 4. **多重检验校正**： - 由于每个基因对应多个变异，需校正假阳性。（即虽有节律性，但并非某一变异带来，而是该基因对应的另一个变异） - 考虑到 **连锁不平衡（LD）**，计算每个基因的 **有效独立变异数（Meff）**。 . Adjusting p values based on Meff可以校正假阳性 - 校正阈值： - 中位数 **Meff = 15**（每个基因）。 0.01/15 约为5 × 10-4，所以 - 设定 **显著性阈值 \(p \leq 5 \times 10^{-4}\)**（相当于Bonferroni校正 \(0.01/15\)）。 - **验证该阈值合理性**： - 对基因型数据进行随机置换（shuffle）随机打乱，发现假阳性率 **< 5%**，支持该阈值的可靠性。 5. **最终筛选标准**： - 保留满足以下条件的变异-基因对（**潜在节律变异-rhyVariant & 节律基因-rhyGene**）： - \(p \leq 5 \times 10^{-4}\) - **峰值-谷值倍数变化（fold change）≥ 1.5**（至少一种基因型）。 --- ### **关键点总结** - **顺式变异范围**：TSS ±1 Mb（涵盖增强子、启动子等调控元件）。 - **节律性分析**：基于 **24小时周期** 的谐波回归，计算振幅和相位。 - **统计严格性**：通过 **Meff校正** 和 **置换检验** 控制假阳性。 - **最终输出**：显著关联的 **rhyQTLs**（影响基因昼夜表达的遗传变异）。 ## 步骤四 节律性变异的评估（图1a步骤4）： 使用dryR包中的drylm函数和贝叶斯信息准则（BIC）评估不同基因型间节律性差异。通过20次降采样和G检验确定差异节律性的显著性，最终定义==rhyVariants的标准==包括：（1）至少在一个基因型亚群中p≤5×10\(^{-4}\)且振幅≥1.5；（2）dryR和G检验显示显著差异节律性；（3）校正后q(BH)<0.05。通过LD区块划分确定独立rhyQTLs。 性能评估：通过随机抽样验证rhyQTLs检测方法的稳健性（补充图4），并利用昼夜RNA-seq数据验证GTEx中鉴定的rhyQTLs（补充图5）。此外，通过NR1D1相关rhyQTLs分析了其对远端基因节律性表达的远程效应（补充图6）。 通路富集分析：使用enrichR包和Enrichr数据库进行，背景基因列表为相应组织中表达的基因（平均计数>10）。显著性阈值设为校正后p<0.05。 GWAS富集分析：从GWAS目录获取SNP-性状/疾病关联，过滤后保留242,822个高质量关联。通过随机匹配SNP数量和MAF分布计算rhyQTLs在GWAS标签SNPs中的富集程度。 LDSC回归：分析15种脂质或脂蛋白相关性状中rhyQTLs和eQTLs的遗传力贡献比例。 共定位分析：使用COLOC包评估rhyQTLs与GWAS信号的共定位，后验概率PP4>0.70视为显著。 VIVA研究GWAS：对815名肥胖易感儿童进行75种肥胖相关性状的GWAS分析，鉴定显著关联（p≤7.7×10\(^{-8}\)）并评估rhyQTLs的富集。 #### 数据与代码可用性 GTEx项目V8的原始基因组测序数据可通过dbGaP获取（phs000424.v8.p2）， 基因表达数据从GTEx门户网站下载。 rhyQTLs的汇总统计信息见补充数据2， 自定义源代码发布于GitHub（https://github.com/YingChen10/rhyQTL）\(^{82}\)。 您的理解基本正确，但需要根据研究目标决定是否需要下载两个网站的数据。以下是具体建议： ### **1. 数据需求分析** | 数据类型 | 存放位置 | 是否需要下载？ | |------------------------|------------------------------|------------------------------------------------------------------------------| | **原始基因组测序数据** | dbGaP (phs000424.v8.p2) | 仅当您需要：<br>• 重新分析原始序列（如call变异）<br>• 验证特定样本的测序质量 | | **基因表达矩阵** | GTEx Portal (bulk_tissue_expression) | 如果只需分析基因表达量、eQTL/rhyQTL等，**只需下载此数据**（已处理好的标准化表达量） | ### **2. 注意事项** #### **（1）GTEx Portal表达数据的优势** - **已预处理**：包含TPM/FPKM标准化表达矩阵，可直接用于分析。 - **样本丰富**：涵盖838名捐赠者的45种组织，附带临床元数据（年龄、性别、死亡时间等）。 - **免审批**：表达数据为开放访问（Open Access），无需申请权限。 #### **（2）dbGaP原始数据的适用场景** - **受限访问**：需提交研究提案并通过dbGaP审批（涉及人类数据隐私）。 - **使用场景举例**：<br> • 开发新的RNA-seq分析流程<br> • 探索GTEx未提供的转录本异构体或罕见变异<br> • 整合其他组学数据（如外显子组测序） --- ### **3. 推荐下载策略** #### **如果您的目标是：** - **分析基因节律性（如本文）** → **仅需GTEx表达数据**（含时间信息文件 [Zenodo链接](https://doi.org/10.5281/zenodo.7215362)）。 - **挖掘新的遗传变异** → 需额外下载dbGaP原始数据（需审批）。 #### **GTEx表达数据下载步骤：** 1. 访问 [GTEx下载页面](https://www.gtexportal.org/home/downloads)。 2. 选择 **"Bulk Tissue Expression"** → 下载 **"Gene TPMs"**（推荐）或 **"Gene Counts"**。 3. 同步下载 **样本属性表**（含组织类型、供体信息等）。 --- ### **4. 数据关联示例** 若您研究rs7740525对AIF1基因的调控： - **GTEx数据**：提供AIF1在各组织中的表达量，验证其是否在心脏中呈现节律性。 - **dbGaP数据**：如需分析该SNP附近的测序覆盖度或寻找共定位变异，则需原始BAM文件。 --- ### **总结** - **常规表达分析**：GTEx Portal数据足够（更快捷）。 - **深度遗传分析**：需结合dbGaP原始数据（流程复杂）。 - 建议先尝试用GTEx数据完成初步分析，必要时再申请dbGaP权限。 #### 参考文献（略） #### 致谢 感谢所有提供宝贵意见和技术支持的同事，以及资助本研究的相关机构。 #### 作者贡献 Y.C.和D.G.设计研究、解读数据并撰写论文，所有作者修订并批准最终稿。Y.C.和P.L.负责基于RNA-seq数据的SNP calling，Y.C.和A.S.处理VIVA项目数据。 #### 竞争利益 作者声明无竞争利益。 #### 附加信息 补充信息及通讯材料请求详见原文链接。 文件名：补充数据1-11 补充数据1. 本研究用于rhyQTL定位的组织样本汇总，包括组织名称及样本量。 补充数据2. 45种组织中rhyQTL的统计摘要。该数据已上传至Zenodo平台（DOI: 10.5281/zenodo.15127913），包含多个按"组织名称.txt"命名的文本文件。每份文件包含19列数据，具体如下： 1. ID - rhy变异体标识符 2. Chromosome - rhyQTL所在染色体 3. Position - rhyQTL位点坐标 4. REF - 参考等位基因 5. ALT - 替代等位基因 6. rhyGene.ID - 节律基因ID 7. rhyGene.name - 节律基因名称 8-10. Sample.size.0/1/2 - 按样本量排序的三个基因型亚组样本量 8. pval_0 - 第一基因型亚组的p值（11-16列为谐波回归拟合参数） 9. phase_0 - 第一亚组的表达峰值相位（24小时制） 10. amp_0 - 第一亚组的对数化峰谷振幅值（log2） 11. pval_1 - 第二基因型亚组的p值 12. phase_1 - 第二亚组的表达峰值相位 13. amp_1 - 第二亚组的对数化峰谷振幅值 14. Chosen.model - DryR选择的模型（取值为2/3/5的整数） 15. G.test_pval - G检验p值 16. pval - 基因型组间节律表达差异的统计学显著性p值 补充数据3. 脑组织中NR1D1相关顺式rhyQTL介导的反式节律基因。 补充数据4. 节律基因峰值相位。每个节律基因的峰值相位定义为表达达到峰值的时间点。若基因存在多个rhyQTL，则取所有对应相位的中位数作为该基因相位。 补充数据5. rhyQTL在基因组调控元件中的富集分析。通过计算观察值与预期值的比值比（odds ratio）评估富集程度，基线区域为人类基因组中所有功能及推定功能区域的联合集。 补充数据6. 45种组织中rhyQTL位点的转录因子结合基序富集分析。通过构建2×2列联表分析转录因子结合基序与遗传变异的共现关系，使用Fisher精确检验计算比值比和p值评估富集显著性。文件按"组织名称.txt"格式存储，包含9列数据： 1. Motif - 基序序列 2. Name - 最佳匹配基序及详细信息 3. TF - 转录因子名称 4-7. 分别表示"rhyQTL+/基序+"、"rhyQTL+/基序-"、"rhyQTL-/基序+"、"rhyQTL-/基序-"的SNP数量 4. OddsRatio - Fisher精确检验计算的比值比 5. p.value - 富集分析显著性p值 补充数据7. rhyQTL相关性状/疾病。列出GWAS中鉴定到≥3个与rhyQTL重叠SNP的性状/疾病。 补充数据8. 人类性状/疾病分类摘要。GWAS Catalog中的性状/疾病按其生理功能分为13大类。 补充数据9. 与15种脂质/脂蛋白相关性状共定位的rhyQTL节律基因。使用COLOC PP4统计方法鉴定与GWAS信号共享变异位点的rhyQTL，并列出影响特定人类表型的节律基因。 补充数据10. VIVA队列中测量的肥胖相关性状。 补充数据11. GWAS鉴定的肥胖相关性状SNP。基于75个肥胖相关表型，按照标准GWAS分析流程（https://github.com/MareesAT/GWA_tutorial）进行分析。 （注：专业术语处理说明： 1. rhyQTL/rhyGene等专业缩写保留不译，通过上下文明确含义 2. DryR/COLOC等分析方法名称保留原名 3. 统计学术语如odds ratio统一译为"比值比" 4. 技术术语如"peak phase"译为"峰值相位"以符合生物学表述习惯 5. 长难句采用拆分法处理，如补充数据6的描述拆分为流程说明和文件结构两部分） File Name: Supplementary Data 1-11 Supplementary Data 1. Summary of tissues used for rhyQTL mapping in this study, including their names and sample sizes. Supplementary Data 2. Summary statistics of rhyQTLs in 45 tissues. The Supplementary Data 2 has been uploaded to Zenodo: https://doi.org/10.5281/zenodo.15127913. The folder contains multiple .txt files named following the pattern tissue_name.txt. Each file includes data organized into 19 columns, which are listed as follows: 1. ID - Identifies the rhyVariant. 2. Chromosome - Indicates the chromosome of rhyQTL. 3. Position - Specifies the position of rhyQTL. 4. REF - Reference allele of rhyQTL. 5. ALT - Alternative allele of rhyQTL. 6. rhyGene.ID - Gene id for the rhyGene. 7. rhyGene.name - Name of the rhyGene. 8-10. Sample.size.0/1/2 - Sample size of the three genotype subgroups, ordered by sample size. 11. pval_0 - p value for the first genotype subgroup. Columns 11-16 are parameters from the harmonic regression fit. 12. phase_0 - Peak phase (24h scale) for the first genotype subgroup. 13. amp_0 - log2 amplitude peak-to-trough for the first genotype subgroup. 14. pval_1 - p value for the second genotype subgroup. 15. phase_1 - Peak phase (24h scale) for the second genotype subgroup. 16. amp_1 - log2 amplitude peak-to trough for the second genotype subgroup. 17. Chosen.model - The selected model given by DryR, with an integer value of 2, 3, or 5. 18. G.test_pval - p value from the G-test. 19. pval - p value reflects the statistical significance of differential rhythmic expression between the two genotype groups. Supplementary Data 3. Trans-rhyGenes mediated by NR1D1-associated cis-rhyQTLs in brain tissues. Supplementary Data 4. Peak phases of rhyGenes. The peak phase of each rhyGene is defined as the time point at which it exhibits the expression peak. For rhyGenes with more than one rhyQTLs, the median value of phases corresponding to all rhyQTLs is taken as the phase of the rhyGene. Supplementary Data 5. Enrichment of rhyQTLs in genomic regulatory elements. The enrichment is calculated as an odds ratio score using the numbers of observed QTLs and expected QTLs located in each annotation category compared to those in baseline regions. Baseline regions are unions of all functional and putative functional regions within the human genome. Supplementary Data 6. TF motif enrichment in genomic loci harboring rhyQTLs across 45 tissues. The relationship between the cooccurrence of transcription factor (TF) binding motifs and genetic variants at different genomic coordinates was analyzed. Genetic variants were categorized based on whether they were identified as rhyQTLs and whether they were located within the motif, forming a 2 × 2 contingency table for each motif. The folder contains multiple .txt files named following the pattern tissue_name.txt. Each file includes data organized into 9 columns, which are listed as follows: 1. Motif: The specific motif sequence being analyzed. 2. Name: The best matching motif and additional details about the motif 3. TF : The name of the transcription factor associated with the motif. 4. rhyQTL+Motif+: The number of SNPs that are both rhyQTLs and located within a motif binding site. 5. rhyQTL+Motif-: The number of SNPs that are rhyQTLs but are not located within a motif binding site. 6. rhyQTL-Motif+: The number of SNPs that are not rhyQTLs but are located within a motif binding site. 7. rhyQTL-Motif-: The number of SNPs that are neither rhyQTLs nor located within a motif binding site. 8. OddsRatio: The odds ratio calculated using Fisher's exact test, indicating the strength of the enrichment. 9. p.value: The p value calculated using Fisher's exact test, indicating the statistical significance of the enrichmentodds ratio was calculated to assess rhyQTL enrichment within motifs, and statistical significance was determined using a two-sided Fisher's exact test. Supplementary Data 7. rhyQTL-associated traits/diseases. Traits/diseases with more than three associated SNPs identified in GWAS that also identified as rhyQTLs are listed. Supplementary Data 8. Summary of human trait or disease categories. Traits/diseases from GWAS Catalog are classified into 13 categories based on their specific physiological functions. Supplementary Data 9. rhyGenes with rhyQTLs that colocalize with loci linked to 15 lipid/lipoprotein-related traits. rhyQTLs sharing variants with significant GWAS signals in 15 lipid/lipoprotein-related traits were identified using COLOC PP4 statistics to estimate their effects on specific human phenotypes. rhyGenes whose rhyQTLs colocalized with GWAS signals for these traits are listed. Supplementary Data 10. Obesity-related traits measured in VIVA cohort. Supplementary Data 11. SNPs associated with obesity-related traits identified by GWAS. GWAS analysis on these 75 obesity-related phenotypes was conducted following the GWAS standard pipeline https://github.com/MareesAT/GWA_tutorial.
